#include "CalculateSunRise.h"



struct SunRiseSetHourMin {
  int hour;
  int minute;
};

//计算日出日落时间
/**
 * (1)先计算出从格林威治时间公元2000年1月1日到计算日天数days;
 * (2)计算从格林威治时间公元2000年1月1日到计算日的世纪数t, 则t＝(days+UTo／360)／36525;
 * (3)计算太阳的平黄径 ： L=280.460+36000.770×t;
 * (4)计算太阳的平近点角 ：G＝357.528+35999.050×t
 * (5)计算太阳的黄道经度 ：λ＝L+1.915×sinG+0.020xsin(2G);
 * (6)计算地球的倾角     ε＝23.4393-0.0130×t;
 * (7)计算太阳的偏差     δ＝arcsin(sinε×sinλ);
 * (8)计算格林威治时间的太阳时间角GHA： GHA=UTo-180-1.915×sinG-0.020×sin(2G)   +2.466×sin(2λ)-0.053×sin(4λ)
 * (9)计算修正值e： e=arcos{[   sinh-sin(Glat)sin(δ)]/cos(Glat)cos(δ)}
 * (10)计算新的日出日落时间 ：UT＝UTo-(GHA+Long±e); 其中“+”表示计算日出时间,“-”表示计算日落时间;
 * (11)比较UTo和UT之差的绝对值,如果大于0.1°即0.007小时,把UT作为新的日出日落时间值,重新从第(2)步开始进行迭代计算,如果UTo和UT之差的绝对值小于0.007小时,则UT即为所求的格林威治日出日落时间;
 * (12)上面的计算以度为单位,即180°=12小时,因此需要转化为以小时表示的时间,再加上所在的时区数Zone,即要计算地的日出日落时间为 ：T=UT/15+Zone
 * 上面的计算日出日落时间方法适用于小于北纬60°和南纬60°之间的区域,如果计算位置为西半球时,经度Long为负数。
 */
// int days_of_month_1[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
// int days_of_month_2[12] = {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
// double h = -0.833;//日出日落时太阳的位置
// double UTo = 180.0;//上次计算的日落日出时间，初始迭代值180.0

// 0、判断是否为闰年：若为闰年，返回1；若不是闰年,返回0
bool CalculateSunRise::leap_year(int year) {
    return (year % 400 == 0) || ((year % 100 != 0) && (year % 4 == 0));
}
// 1、求从格林威治时间公元2000年1月1日到计算日天数days
int CalculateSunRise::days(int year, int month, int date) {
  int days_of_month_1[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  int days_of_month_2[12] = {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
    int i, a = 0;
    for (i = 2000; i < year; i++) {
        if (leap_year(i)) a = a + 366;
        else a = a + 365;
    }

    if (leap_year(year)) {
        for (i = 0; i < month - 1; i++) {
            a = a + days_of_month_2[i];
        }
    } else {
        for (i = 0; i < month - 1; i++) {
            a = a + days_of_month_1[i];
        }
    }
    a = a + date;
    return a;
}
//求格林威治时间公元2000年1月1日到计算日的世纪数t
double CalculateSunRise::t_century(int days, double uto) {
    return (double)(days + uto / 360) / 36525;
}
//求太阳的平黄径
double CalculateSunRise::L_sun(double t_century) {
    return (280.460 + 36000.770 * t_century);
}
//求太阳的平近点角
double CalculateSunRise::G_sun(double t_century) {
    return (357.528 + 35999.050 * t_century);
}
//求黄道经度
double CalculateSunRise::ecliptic_longitude(double L_sun, double G_sun) {
    return (L_sun + 1.915 * sin(G_sun * PI / 180) + 0.02 * sin(2 * G_sun * PI / 180));
}
//求地球倾角
double CalculateSunRise::earth_tilt(double t_century) {
    return (23.4393 - 0.0130 * t_century);
}
//求太阳偏差
double CalculateSunRise::sun_deviation(double earth_tilt, double ecliptic_longitude) {
    return (180 / PI * asin(sin(PI / 180 * earth_tilt) * sin(PI / 180 * ecliptic_longitude)));
}
//求格林威治时间的太阳时间角GHA
double CalculateSunRise::GHA(double UTo, double G_sun, double ecliptic_longitude) {
    return (UTo - 180 - 1.915 * sin(G_sun * PI / 180) - 0.02 * sin(2 * G_sun * PI / 180)
            + 2.466 * sin(2 * ecliptic_longitude * PI / 180)
            - 0.053 * sin(4 * ecliptic_longitude * PI / 180));
}
//求修正值e
double CalculateSunRise::e(double h, double glat, double sun_deviation) {
    return 180 / PI * acos((sin(h * PI / 180) - sin(glat * PI / 180)
                                                                    * sin(sun_deviation * PI / 180)) / (cos(glat * PI / 180)
                                                                                                                  * cos(sun_deviation * PI / 180)));
}
//求日出时间
double CalculateSunRise::UT_rise(double UTo, double GHA, double glong, double e) {
    return (UTo - (GHA + glong + e));
}
//求日落时间
double CalculateSunRise::UT_set(double UTo, double GHA, double glong, double e) {
    return (UTo - (GHA + glong - e));
}
//判断并返回结果（日出）
double CalculateSunRise::result_rise(double UT, double UTo, double glong, double glat, int year, int month, int date) {
    double d;
    if (UT >= UTo) d = UT - UTo;
    else d = UTo - UT;
    if (d >= 0.1) {
        UTo = UT;
        UT = UT_rise(UTo,
                     GHA(UTo, G_sun(t_century(days(year, month, date), UTo)),
                         ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                            G_sun(t_century(days(year, month, date), UTo)))),
                     glong,
                     e(h, glat, sun_deviation(earth_tilt(t_century(days(year, month, date), UTo)),
                                              ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                                                 G_sun(t_century(days(year, month, date), UTo))))));
        result_rise(UT, UTo, glong, glat, year, month, date);
    }

    return UT;

}

//判断并返回结果（日落）
double CalculateSunRise::result_set(double UT, double UTo, double glong, double glat, int year, int month, int date) {
    double d;
    if (UT >= UTo) d = UT - UTo;
    else d = UTo - UT;
    if (d >= 0.1) {
        UTo = UT;
        UT = UT_set(UTo,
                    GHA(UTo, G_sun(t_century(days(year, month, date), UTo)),
                        ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                           G_sun(t_century(days(year, month, date), UTo)))),
                    glong,
                    e(h, glat, sun_deviation(earth_tilt(t_century(days(year, month, date), UTo)),
                                             ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                                                G_sun(t_century(days(year, month, date), UTo))))));
        result_set(UT, UTo, glong, glat, year, month, date);
    }
    return UT;
}

// 日出
SunRiseSetHourMin CalculateSunRise::getSunrise(double longitude, double latitude, int year, int month, int date) {
   double sunrise = result_rise(UT_rise(UTo,
                                  GHA(UTo, G_sun(t_century(days(year, month, date), UTo)),
                                      ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                                         G_sun(t_century(days(year, month, date), UTo)))),
                                  longitude,
                                  e(h, latitude, sun_deviation(earth_tilt(t_century(days(year, month, date), UTo)),
                                                               ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                                                                  G_sun(t_century(days(year, month, date), UTo)))))), UTo, longitude, latitude, year, month, date);
    // return (int)(sunrise / 15 + 8) + ":" + (int) (60 * (sunrise / 15 + 8 - (int) (sunrise / 15 + 8)));
   struct SunRiseSetHourMin sunRiseHourMin; 
   sunRiseHourMin.hour = (int)(sunrise / 15 + 8);
   sunRiseHourMin.minute = (int) (60 * (sunrise / 15 + 8 - (int) (sunrise / 15 + 8)));
   return sunRiseHourMin;
}

// 日落
SunRiseSetHourMin CalculateSunRise::getSunset(double longitude, double latitude, int year, int month, int date) {
    double sunset = result_set(UT_set(UTo,
                               GHA(UTo, G_sun(t_century(days(year, month, date), UTo)),
                                   ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                                      G_sun(t_century(days(year, month, date), UTo)))),
                               longitude,
                               e(h, latitude, sun_deviation(earth_tilt(t_century(days(year, month, date), UTo)),
                                                            ecliptic_longitude(L_sun(t_century(days(year, month, date), UTo)),
                                                                               G_sun(t_century(days(year, month, date), UTo)))))), UTo, longitude, latitude, year, month, date);
    // return (int) (sunset / 15 + 8) + ":" + (int) (60 * (sunset / 15 + 8 - (int) (sunset / 15 + 8)));
    struct SunRiseSetHourMin sunSetHourMin;
    sunSetHourMin.hour = (int) (sunset / 15 + 8);
    sunSetHourMin.minute = (int) (60 * (sunset / 15 + 8 - (int) (sunset / 15 + 8)));
    return sunSetHourMin;
}

//判断白天
bool CalculateSunRise::bDaytime(double lon, double lat, int year, int month, int date, int hour, int minute){
    initClass();
    //Serial.println(F("lon lat year month date hour minute: ") + (String)lon + F(",") + lat + F(",") + year + F(",") + month + F(",") + date + F(",") + hour + F(",") + minute);
    struct SunRiseSetHourMin sunsetHourMin = getSunset(lon, lat, year,month,date);
    struct SunRiseSetHourMin sunriseHourMin = getSunrise(lon,lat,year,month,date);
    long bDaytime_tmpS = 3600 * (long)hour + 60 * (long)minute;
    long bDaytime_sunrise = 3600 * (long)sunriseHourMin.hour + 60 * (long)sunriseHourMin.minute;
    long bDaytime_sunset = 3600 * (long)sunsetHourMin.hour + 60 * (long)sunsetHourMin.minute;
    return bDaytime_tmpS > bDaytime_sunrise && bDaytime_tmpS < bDaytime_sunset;
}
void CalculateSunRise::initClass() {
  h = -0.833;//日出日落时太阳的位置
  UTo = 180.0;//上次计算的日落日出时间，初始迭代值180.0
}